function stamstat()

bands={'theta','alpha','beta','gamma'};

inDir='C:\Users\matthew\Desktop\Projects\Batch\analysis\mc\ascii\stam\';
outDir='C:\Users\matthew\Desktop\Projects\Batch\analysis\mc\ascii\stam\';


mlc=1:15;
mlf=16:31;
mlo=32:41;
mlp=42:50;
mlt=51:71;

mrc=72:86;
mrf=87:102;
mro=103:112;
mrp=113:121;
mrt=122:142;

labels={'mlc','mlf','mlo','mlp','mlt','mrc','mrf','mro','mrp','mrt'};
areas={mlc;mlf;mlo;mlp;mlt;mrc;mrf;mro;mrp;mrt};
alphalevel=0.05;

filter='plistam*'
files=dir(fullfile(inDir,filter));

for j=1:length(files)  
    load(fullfile(inDir,files(j).name),'subj')
    nSubjs=size(subj.one,1);
    for i=1:nSubjs

            statsingle=statOverTr(squeeze(subj.one(i,:,:,:)),squeeze(subj.two(i,:,:,:)),areas,alphalevel);
            name=subj.oneLabels{i};    
            [r,subjName]=strtok(name{1},'-');
            [band,rr]=strtok(r,'_');
            statfile=strcat('stat_',band);
            statfile=strcat(statfile,subjName);

            save(fullfile(outDir,statfile),'statsingle')
    end
end

for i=1:length(bands)
    statOverSubjs(outDir,outDir,bands{i},alphalevel)
end



function [res]=statOverTr(one_pliXT,two_pliXT,areas,alphalevel)

nTrial=min(size(one_pliXT,1),size(two_pliXT,1));
one_trialXintra_area=zeros(nTrial,length(areas));
two_trialXintra_area=zeros(nTrial,length(areas));

interAreaLen=((length(areas)*length(areas))-length(areas))/2;

one_trialXinter_area=zeros(nTrial,interAreaLen);
two_trialXinter_area=zeros(nTrial,interAreaLen);


s=size(one_pliXT);
indx=ones(s(2),s(3));
indx=indx-tril(indx);
indx=logical(indx);

one_avg_pliXtrial=zeros(nTrial,1);
two_avg_pliXtrial=zeros(nTrial,1);


for i=1:nTrial
   
   tone=squeeze(one_pliXT(i,:,:));
   ttwo= squeeze(two_pliXT(i,:,:));
    
   one_trialXintra_area(i,:)=shortPLI(tone,areas);
   two_trialXintra_area(i,:)=shortPLI(ttwo,areas);
   
   one_trialXinter_area(i,:)=longPLI(tone,areas);
   two_trialXinter_area(i,:)=longPLI(ttwo,areas);
   
   one_avg_pliXtrial(i)=mean(tone(indx));
   two_avg_pliXtrial(i)=mean(ttwo(indx));

end

res.intra={one_trialXintra_area,two_trialXintra_area};
res.inter={one_trialXinter_area,two_trialXinter_area};
res.avg={one_avg_pliXtrial,two_avg_pliXtrial};

[res.h_intra,res.p_intra]=ttest(one_trialXintra_area,two_trialXintra_area,alphalevel);
[res.h_inter,res.p_inter]=ttest(one_trialXinter_area,two_trialXinter_area,alphalevel);
[res.h_avgpli,res.p_avgpli]=ttest(one_avg_pliXtrial,two_avg_pliXtrial,alphalevel);






function statOverSubjs(inDir,outDir,bandName,alphalevel)

filter=strcat('stat_*',bandName);
filter=strcat(filter,'*');

files=dir(fullfile(inDir,filter));
if(length(files)==0)return;end 
nSubjs=length(files);

load(fullfile(inDir,files(1).name),'statsingle');

nAreas=size(statsingle.intra{1},2)
nPairAreas=size(statsingle.inter{1},2);

one_SubjXintraA=zeros(nSubjs,nAreas);
two_SubjXintraA=zeros(nSubjs,nAreas);
one_SubjXinterA=zeros(nSubjs,nPairAreas);
two_SubjXinterA=zeros(nSubjs,nPairAreas);
one_SubjAVG=zeros(nSubjs,1);
two_SubjAVG=zeros(nSubjs,1);

for i=1:nSubjs
    load(fullfile(inDir,files(i).name),'statsingle');

    one_SubjXintraA(i,:)=mean(statsingle.intra{1},1);
    two_SubjXintraA(i,:)=mean(statsingle.intra{2},1);
    one_SubjXinterA(i,:)=mean(statsingle.inter{1},1);
    two_SubjXinterA(i,:)=mean(statsingle.inter{2},1);
    one_SubjAVG(i)=mean(statsingle.avg{1},1);
    two_SubjAVG(i)=mean(statsingle.avg{2},1);
    
    [r,subjName]=strtok(files(i).name,'-');
    subjName=subjName(2:end);
    subjName=strtok(subjName,'.');
    multi.subjNames{i}=subjName;  
end

multi.intra={one_SubjXintraA,two_SubjXintraA};
multi.inter={one_SubjXinterA,two_SubjXinterA};
multi.avg={one_SubjAVG,two_SubjAVG};

[multi.h_intra,multi.p_intra]=ttest(one_SubjXintraA,two_SubjXintraA,alphalevel);
[multi.h_inter,multi.p_inter]=ttest(one_SubjXinterA,two_SubjXinterA,alphalevel);
[multi.h_avgPLI,multi.p_avgPLI]=ttest(one_SubjAVG,two_SubjAVG,alphalevel);

outfile=fullfile(outDir,strcat('statMultiSubjs_',bandName));

save(outfile,'multi');





function [v]=shortPLI(pli,areas)

v=zeros(1,length(areas));
for i=1:length(areas)
    m=pli(areas{i},areas{i}); 
    indx=ones(size(v));
    indx=indx-tril(indx);
    indx=logical(indx);    
    v(i)=mean(m(indx));
end

function [interArea]=longPLI(pli,areas)

v=zeros(length(areas),length(areas));
for i=1:length(areas)-1
    for j=i+1:length(areas)
        m=pli(areas{i},areas{j});
        s=size(m);
        v(i,j)=mean(reshape(m,s(1)*s(2),1));
    end
end

indx=ones(size(v));
indx=indx-tril(indx);
indx=logical(indx);
%indx=logical(indx);%%warning if phase is 0
interArea=v(indx)';








